home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / ureadseq.c < prev    next >
Text File  |  1996-07-05  |  65KB  |  2,303 lines
  1. /* File: ureadseq.c
  2.  *
  3.  * Reads and writes nucleic/protein sequence in various
  4.  * formats. Data files may have multiple sequences.
  5.  *
  6.  * Copyright 1990 by d.g.gilbert
  7.  * biology dept., indiana university, bloomington, in 47405
  8.  * e-mail: gilbertd@bio.indiana.edu
  9.  *
  10.  * This program may be freely copied and used by anyone.
  11.  * Developers are encourged to incorporate parts in their
  12.  * programs, rather than devise their own private sequence
  13.  * format.
  14.  *
  15.  * This should compile and run with any ANSI C compiler.
  16.  *
  17.  */
  18.  
  19.  
  20. #include <stdio.h>
  21. #include <ctype.h>
  22. #include <string.h>
  23. #include <stdlib.h>
  24.  
  25. #define UREADSEQ_G
  26. #include "ureadseq.h"
  27.  
  28.  
  29.  
  30. #ifdef DCLAP
  31. #include <ncbi.h> /* for ncbimem.h */
  32.  
  33. #    define MALLOC(size)            (char*)Nlm_MemNew(size)
  34. #    define REALLOC(p,size)     p= (char*)Nlm_MemMore(p,size)
  35. # define FREE(p)                    Nlm_MemFree(p)
  36. # define MEMCPY(d,t,n)        Nlm_MemCopy( d, t, n)
  37.  
  38. #else
  39.  
  40. #    define MALLOC(size)            (char*)malloc(size)
  41. #    define REALLOC(p,size)     p= (char*)realloc(p,size)
  42. # define FREE(p)                    free(p)
  43. # define MEMCPY(d,t,n)        memcpy( d, t, n)
  44. #endif
  45.  
  46. #ifndef Local
  47. # define Local      static    /* local functions */
  48. #endif
  49.  
  50. #define kStartLength 10240     /* 500  500000 to temp fix Unix bug */
  51. long     aStartLength = kStartLength;
  52.  
  53. const char *aminos      = "ABCDEFGHIKLMNPQRSTVWXYZ*";
  54. const char *primenuc    = "ACGTU";
  55. const char *protonly    = "EFIPQZ";
  56.  
  57. const char kNocountsymbols[5]  = "_.-?";
  58. const char stdsymbols[6]  = "_.-*?";
  59. const char allsymbols[32] = "_.-*?<>{}[]()!@#$%^&=+;:'/|`~\"\\";
  60. static const char *seqsymbols   = allsymbols;
  61.  
  62. const char nummask[11]   = "0123456789";
  63. const char nonummask[11] = "~!@#$%^&*(";
  64.  
  65. FILE    * gFile;
  66. long    gLinestart = 0;
  67.  
  68.  
  69.  
  70. /*
  71.     use general form of isseqchar -- all chars + symbols.
  72.     no formats except nbrf (?) use symbols in data area as
  73.     anything other than sequence chars.
  74. */
  75.  
  76.  
  77.  
  78.                           /* Local variables for readSeq: */
  79. struct ReadSeqVars {
  80.   short choice, err, nseq, memstep;
  81.   long  seqlen, maxseq, seqlencount, estlen;
  82.   short topnseq;
  83.   long  topseqlen;
  84.   const char *fname;
  85.   char *seq, *seqid, matchchar;
  86.   boolean allDone, done, filestart, addit;
  87.   FILE  *f;
  88.   long  linestart;
  89.   char  s[256], *sp;
  90.  
  91.   int (*isseqchar)(int);  
  92.     ReadWriteProc    reader;
  93. };
  94.  
  95.  
  96.  
  97. int Strcasecmp(const char *a, const char *b)  /* from Nlm_StrICmp */
  98. {
  99.   int diff, done;
  100.   if (a == b)  return 0;
  101.   done = 0;
  102.   while (! done) {
  103.     diff = to_upper(*a) - to_upper(*b);
  104.     if (diff) return diff;
  105.     if (*a == '\0') done = 1;
  106.     else { a++; b++; }
  107.     }
  108.   return 0;
  109. }
  110.  
  111. int Strncasecmp(const char *a, const char *b, long maxn) /* from Nlm_StrNICmp */
  112. {
  113.   int diff, done;
  114.   if (a == b)  return 0;
  115.   done = 0;
  116.   while (! done) {
  117.     diff = to_upper(*a) - to_upper(*b);
  118.     if (diff) return diff;
  119.     if (*a == '\0') done = 1;
  120.     else {
  121.       a++; b++; maxn--;
  122.       if (! maxn) done = 1;
  123.       }
  124.     }
  125.   return 0;
  126. }
  127.  
  128.  
  129. #define EOFreader(V)    (boolean)(*V->reader)((char *)NULL, 0L, kRSFile_Eof)
  130. #define REWINDreader(V)    (*V->reader)((char *)NULL, 0L, kRSFile_Rewind)
  131.  
  132. Local long  readWriteFromFile(char* line, long maxline, short action)
  133. {
  134.     static boolean geof = false;
  135.     long   bytes;
  136.     /*fprintf(stderr," r/w action= %d,  maxline= %ld ", action, maxline);*/
  137.     /* ! one case where _Read returns 0 bytes at end of file, but _Eof returns false !! */
  138.     
  139.     switch (action) {
  140.     
  141.     case kRSFile_Eof:
  142.         if (geof) {  /* fix for bad feof() ! */
  143.             char ch= fgetc( gFile);
  144.             if (ch != EOF) { ungetc( ch, gFile); geof= false; }
  145.             }
  146.         else geof= feof( gFile);
  147.         return geof;
  148.  
  149.     case kRSFile_Read:
  150.         bytes= (long) fgets (line, (int)maxline, gFile);
  151.         if (!bytes) geof= true; /* fix for bad feof() ! */
  152.         return bytes;
  153.     
  154.     case kRSFile_Write:
  155.         geof= false;
  156.         fputs( line, gFile);
  157.         return 0;
  158.  
  159.     case kRSFile_Seek:
  160.         geof= false;
  161.         fseek( gFile, maxline, 0);
  162.         return 0;
  163.  
  164.     case kRSFile_End:
  165.         fseek( gFile, 0, 2);
  166.         geof= true;
  167.         return 0;
  168.  
  169.     case kRSFile_Rewind:
  170.         rewind(gFile);
  171.         geof= false;
  172.         return 0;
  173.  
  174.     case kRSFile_Tell:
  175.         return ftell(gFile);
  176.         
  177.     default:
  178.         return 0;
  179.     }
  180. }
  181.  
  182.  
  183.  
  184. enum seqflag { 
  185.     primenucflag = 1,
  186.     dnanucflag = 2,
  187.     rnanucflag = 4,
  188.     aminoflag = 8,
  189.     protonlyflag = 0x10,
  190.     seqsymflag = 0x20
  191.     };
  192.  
  193. char seqflags[128] = {
  194.   0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,
  195.   0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,
  196. /*    ,  ! ,  " ,  # ,  $ ,  % ,  & ,  ' ,  ( ,  ) ,  * ,  + ,  , ,  - ,  . ,  / ,*/
  197.   0x0 ,0x20,0x20,0x20,0x0 ,0x0 ,0x20,0x20,0x20,0x20,0x28,0x20,0x0 ,0x20,0x20,0x20,
  198. /*  0 ,  1 ,  2 ,  3 ,  4 ,  5 ,  6 ,  7 ,  8 ,  9 ,  : ,  ; ,  < ,  = ,  > ,  ? ,*/
  199.   0x0, 0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x0 ,0x20,0x20,0x20,0x20,0x20,0x20,
  200. /*  @ ,  A ,  B ,  C ,  D ,  E ,  F ,  G ,  H ,  I ,  J ,  K ,  L ,  M ,  N ,  O ,*/
  201.   0x20,0x9 ,0x8 ,0x9 ,0x8 ,0x18,0x18,0x9 ,0x8 ,0x18,0x0 ,0x8 ,0x8 ,0x8 ,0x8 ,0x0 ,
  202. /*  P ,  Q ,  R ,  S ,  T ,  U ,  V ,  W ,  X ,  Y ,  Z ,  [ ,  \ ,  ] ,  ^ ,  _ ,*/
  203.   0x18,0x18,0x8 ,0x8 ,0xB ,0x5 ,0x8 ,0x8 ,0x8 ,0x8 ,0x18,0x20,0x20,0x20,0x20,0x20,
  204. /*  ` ,  a ,  b ,  c ,  d ,  e ,  f ,  g ,  h ,  i ,  j ,  k ,  l ,  m ,  n ,  o ,*/
  205.   0x20,0x9 ,0x8 ,0x9 ,0x8 ,0x18,0x18,0x9 ,0x8 ,0x18,0x0 ,0x8 ,0x8 ,0x8 ,0x8 ,0x0 ,
  206. /*  p ,  q ,  r ,  s ,  t ,  u ,  v ,  w ,  x ,  y ,  z ,  { ,  | ,  } ,  ~ ,   ,*/
  207.   0x18,0x18,0x8 ,0x8 ,0xB ,0x5 ,0x8 ,0x8 ,0x8 ,0x8 ,0x18,0x20,0x20,0x20,0x20,0x0 
  208. };
  209.  
  210. boolean isprimenuc(char c) { return (c>127||c<0)?false:seqflags[c] & primenucflag; }
  211. boolean isdnanuc(char c) { return (c>127||c<0)?false:seqflags[c] & dnanucflag; }
  212. boolean isrnanuc(char c) { return (c>127||c<0)?false:seqflags[c] & rnanucflag; }
  213. boolean isamino(char c) { return (c>127||c<0)?false:seqflags[c] & aminoflag; }
  214. boolean isprotonly(char c) { return (c>127||c<0)?false:seqflags[c] & protonlyflag; }
  215. boolean isseqsym(char c) { return (c>127||c<0)?false:seqflags[c] & seqsymflag; }
  216.  
  217.  
  218. int isSeqChar(int c)
  219. {
  220.  /*  return (isalpha(c) || strchr(seqsymbols,c)); */
  221.   if (isalpha(c) || isseqsym(c)) return c;
  222.   else return 0;
  223. }
  224.  
  225. int isGCGSeqChar(int c)
  226. {
  227.  /*  return (isalpha(c) || strchr(seqsymbols,c)); */
  228.   if (isalpha(c) || isseqsym(c))  {
  229.     if (c == '.') return '-'; /* do the indel translate */
  230.     else return c;
  231.     }
  232.   else return 0;
  233. }
  234.  
  235. int isSeqNumChar(int c)
  236. {
  237.   /* return (isalnum(c) || strchr(seqsymbols,c)); */
  238.   if (isalnum(c) || isseqsym(c)) return c;
  239.   else return 0;
  240. }
  241.  
  242.  
  243. int isAnyChar(int c)
  244. {
  245.  /*  return isprint(c); /* wrap in case isascii is macro */
  246.  if (isprint(c)) return c;
  247.  else return 0;
  248. }
  249.  
  250.  
  251. Local void callreadline( ReadWriteProc reader, char *s)
  252. {
  253.   char  *cp;
  254.  
  255.     /*fprintf(stderr,"callreadline &reader= %ld \n", reader);*/
  256.     
  257.     gLinestart= (*reader)( (char*)NULL, 0L, kRSFile_Tell);
  258.     if (0 == (*reader)( s, 256L, kRSFile_Read) )
  259.         *s = 0;
  260.   else {
  261.     cp = strchr(s, '\n');
  262.     if (cp != NULL) *cp = 0;
  263.     }
  264. }
  265.  
  266. /********
  267. Local void readline(FILE *f, char *s, long *linestart)
  268. {
  269.   char  *cp;
  270.  
  271.     *linestart= ftell(f);
  272.     if (NULL == fgets(s, 256, f))
  273.     *s = 0;
  274.   else {
  275.     cp = strchr(s, '\n');
  276.     if (cp != NULL) *cp = 0;
  277.     }
  278. }
  279. ******/
  280.  
  281. Local void getline(struct ReadSeqVars *V)
  282. {
  283.   /*readline(V->f, V->s, &V->linestart);  */
  284.     callreadline(V->reader, V->s);
  285. }
  286.  
  287. Local void ungetline(struct ReadSeqVars *V)
  288. {
  289.     /*fseek(V->f, V->linestart, 0); */
  290.     (void) (*(V->reader))( (char*)NULL, gLinestart, kRSFile_Seek);  
  291. }
  292.  
  293.  
  294. Local void addseq(char *s, struct ReadSeqVars *V)
  295. {
  296.   char  *ptr;
  297.   register char seqc;
  298.   
  299.   if (V->addit) while (*s != 0) {
  300.     if ((seqc= (V->isseqchar)(*s)) != 0) {
  301.       if (V->seqlen >= V->maxseq) {
  302.         V->maxseq += aStartLength;
  303.         V->memstep++;
  304. /* revise this to 
  305.   (a) store all of long seq on disk file (tmpnam()),
  306.         -- offload to disk instead of realloc until done reading
  307.   (b) increase aStartLength if seq is getting large -- larger chunksize
  308.   (c) check ftell(V->f) to get max seq size, reduce if see it is multiseq file?
  309. */
  310. #if 1
  311.                 if (V->memstep % 10 == 0) {
  312.                     /* fix for memory fragmenting... write to disk, free, then malloc. */
  313.                     FILE * tmpf= tmpfile();
  314.                     fwrite( V->seq, 1, V->seqlen+1, tmpf);
  315.                     FREE(V->seq);
  316.                     ptr = MALLOC(V->maxseq+1);
  317.                     rewind( tmpf);
  318.                     fread( ptr, 1, V->seqlen+1, tmpf);
  319.                     fclose( tmpf); /* deleted here !? */
  320.                     V->seq= ptr;
  321.                     }
  322.                 else {
  323.                     ptr = V->seq;
  324.                 REALLOC(ptr, V->maxseq+1); 
  325.                 }
  326. #else
  327.                 ptr = V->seq;
  328.             REALLOC(ptr, V->maxseq+1); /* mac fills mem - not releasing old ptr !? */
  329. #endif
  330.         if (ptr==NULL) {
  331.           V->err = eMemFull;
  332.           return;
  333.           }
  334.         else V->seq = ptr;
  335.                 }
  336. #ifdef TRANSLATE
  337.       V->seq[(V->seqlen)++] = gTranslate[seqc];
  338. #else
  339.       V->seq[(V->seqlen)++] = seqc;
  340. #endif
  341.       }
  342.     s++;
  343.     }
  344. }
  345.  
  346. Local void countseq(char *s, struct ReadSeqVars *V)
  347.  /* this must count all valid seqq chars, for some formats (paup-sequential) even
  348.     if we are skipping seqq... */
  349. {
  350.   while (*s != 0) {
  351.     if ((V->isseqchar)(*s)) {
  352.       (V->seqlencount)++;
  353.       }
  354.     s++;
  355.     }
  356. }
  357.  
  358.  
  359. Local void addinfo(char *s, struct ReadSeqVars *V)
  360. {
  361.   char s2[256], *si;
  362.   boolean saveadd;
  363.   int (*oldIsSeqChar)(int);  
  364.   
  365.   si = s2;
  366.   while (*s == ' ') s++;
  367.   sprintf(si, " %d)  %s\n", V->nseq, s);
  368.  
  369.   saveadd = V->addit;
  370.   V->addit = true;
  371.   oldIsSeqChar= V->isseqchar;
  372.   V->isseqchar = isAnyChar;
  373.   addseq( si, V);
  374.   V->addit = saveadd;
  375.   V->isseqchar = oldIsSeqChar;  
  376. }
  377.  
  378.  
  379.  
  380.  
  381. Local void readLoop(short margin, boolean addfirst,
  382.             boolean (*endTest)(boolean *addend, boolean *ungetend, struct ReadSeqVars *V),
  383.             struct ReadSeqVars *V)
  384. {
  385.   boolean addend = false;
  386.   boolean ungetend = false;
  387.  
  388.   V->nseq++;
  389.   if (V->choice == kListSequences) V->addit = false;
  390.   else V->addit = (V->nseq == V->choice);
  391.   if (V->addit) V->seqlen = 0;
  392.  
  393.   if (addfirst) addseq(V->s, V);
  394.   do {
  395.     getline(V);
  396.     V->done = EOFreader(V); /* feof(V->f); */
  397.     V->done |= (*endTest)( &addend, &ungetend, V);
  398.     if (V->addit && (addend || !V->done) && (strlen(V->s) > margin)) {
  399.       addseq( (V->s)+margin, V);
  400.     }
  401.   } while (!V->done);
  402.  
  403.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  404.   else {
  405.     V->allDone = (V->nseq >= V->choice);
  406.     if (V->allDone && ungetend) ungetline(V);
  407.     }
  408. }
  409.  
  410.  
  411.  
  412. Local boolean endIG( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  413. {
  414.   *addend = true; /* 1 or 2 occur in line w/ bases */
  415.   *ungetend= false;
  416.   return((strchr(V->s,'1')!=NULL) || (strchr(V->s,'2')!=NULL));
  417. }
  418.  
  419. Local void readIG(struct ReadSeqVars *V)
  420. {
  421. /* 18Aug92: new IG format -- ^L between sequences in place of ";" */
  422.   char  *si;
  423.  
  424.   while (!V->allDone) {
  425.     do {
  426.       getline(V);
  427.       for (si= V->s; *si != 0 && *si < ' '; si++) *si= ' '; /* drop controls */
  428.       if (*si == 0) *V->s= 0; /* chop line to empty */
  429.     } while (! (EOFreader(V) || ((*V->s != 0) && (*V->s != ';') ) ));
  430.     if (EOFreader(V))
  431.       V->allDone = true;
  432.     else {
  433.       strcpy(V->seqid, V->s);
  434.       readLoop(0, false, endIG, V);
  435.       }
  436.   }
  437. }
  438.  
  439.  
  440.  
  441. Local boolean endStrider( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  442. {
  443.   *addend = false;
  444.   *ungetend= false;
  445.   return (strstr( V->s, "//") != NULL);
  446. }
  447.  
  448. Local void readStrider(struct ReadSeqVars *V)
  449. { /* ? only 1 seqq/file ? */
  450.  
  451.   while (!V->allDone) {
  452.     getline(V);
  453.     if (strstr(V->s,"; DNA sequence  ") == V->s)
  454.       strcpy(V->seqid, (V->s)+16);
  455.     else
  456.       strcpy(V->seqid, (V->s)+1);
  457.     while ((!EOFreader(V)) && (*V->s == ';')) {
  458.       getline(V);
  459.       }
  460.     if (EOFreader(V)) V->allDone = true;
  461.     else readLoop(0, true, endStrider, V);
  462.   }
  463. }
  464.  
  465.  
  466. Local boolean endPIR( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  467. {
  468.   *addend = false;
  469.   *ungetend= (strstr(V->s,"ENTRY") == V->s);
  470.   return ((strstr(V->s,"///") != NULL) || *ungetend);
  471. }
  472.  
  473. Local void readPIR(struct ReadSeqVars *V)
  474. { /*PIR -- many seqs/file */
  475.  
  476.   while (!V->allDone) {
  477.     while (! (EOFreader(V) 
  478.                 || strstr(V->s,"ENTRY")  || strstr(V->s,"SEQUENCE")) )
  479.       getline(V);
  480.     strcpy(V->seqid, (V->s)+16);
  481.     while (! (EOFreader(V) || strstr(V->s,"SEQUENCE") == V->s))
  482.       getline(V);
  483.     readLoop(0, false, endPIR, V);
  484.  
  485.     if (!V->allDone) {
  486.      while (! (EOFreader(V) || ((*V->s != 0)
  487.        && (strstr( V->s,"ENTRY") == V->s))))
  488.         getline(V);
  489.       }
  490.     if (EOFreader(V)) V->allDone = true;
  491.   }
  492. }
  493.  
  494.  
  495. Local boolean endGB( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  496. {
  497.   *addend = false;
  498.   *ungetend= (strstr(V->s,"LOCUS") == V->s);
  499.   return ((strstr(V->s,"//") != NULL) || *ungetend);
  500. }
  501.  
  502. Local void readGenBank(struct ReadSeqVars *V)
  503. { /*GenBank -- many seqs/file */
  504.  
  505.   while (!V->allDone) {
  506.     strcpy(V->seqid, (V->s)+12);
  507.     while (! (EOFreader(V) || strstr(V->s,"ORIGIN") == V->s))
  508.       getline(V);
  509.     readLoop(0, false, endGB, V);
  510.  
  511.     if (!V->allDone) {
  512.      while (! (EOFreader(V) || ((*V->s != 0)
  513.        && (strstr( V->s,"LOCUS") == V->s))))
  514.         getline(V);
  515.       }
  516.     if (EOFreader(V)) V->allDone = true;
  517.   }
  518. }
  519.  
  520.  
  521. Local boolean endNBRF( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  522. {
  523.   char  *a;
  524.  
  525.   if ((a = strchr(V->s, '*')) != NULL) { /* end of 1st seqq */
  526.     /* "*" can be valid base symbol, drop it here */
  527.     *a = 0;
  528.     *addend = true;
  529.     *ungetend= false;
  530.     return(true);
  531.     }
  532.   else if (*V->s == '>') { /* start of next seqq */
  533.     *addend = false;
  534.     *ungetend= true;
  535.     return(true);
  536.     }
  537.   else
  538.     return(false);
  539. }
  540.  
  541. Local void readNBRF(struct ReadSeqVars *V)
  542. {
  543.   while (!V->allDone) {
  544.     strcpy(V->seqid, (V->s)+4);
  545.     getline(V);   /*skip title-junk line*/
  546.     readLoop(0, false, endNBRF, V);
  547.     if (!V->allDone) {
  548.      while (!(EOFreader(V) || (*V->s != 0 && *V->s == '>')))
  549.         getline(V);
  550.       }
  551.     if (EOFreader(V)) V->allDone = true;
  552.   }
  553. }
  554.  
  555.  
  556.  
  557. Local boolean endPearson( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  558. {
  559.   *addend = false;
  560.   *ungetend= true;
  561.   return(*V->s == '>');
  562. }
  563.  
  564. Local void readPearson(struct ReadSeqVars *V)
  565. {
  566.   while (!V->allDone) {
  567.     strcpy(V->seqid, (V->s)+1);
  568.     readLoop(0, false, endPearson, V);
  569.     if (!V->allDone) {
  570.      while (!(EOFreader(V) || ((*V->s != 0) && (*V->s == '>'))))
  571.         getline(V);
  572.       }
  573.     if (EOFreader(V)) V->allDone = true;
  574.   }
  575. }
  576.  
  577.  
  578.  
  579. Local boolean endEMBL( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  580. {
  581.   *addend = false;
  582.   *ungetend= (strstr(V->s,"ID   ") == V->s);
  583.   return ((strstr(V->s,"//") != NULL) || *ungetend);
  584. }
  585.  
  586. Local void readEMBL(struct ReadSeqVars *V)
  587. {
  588.   while (!V->allDone) {
  589.     strcpy(V->seqid, (V->s)+5);
  590.     do {
  591.       getline(V);
  592.     } while (!(EOFreader(V) || (strstr(V->s,"SQ   ") == V->s)));
  593.  
  594.     readLoop(0, false, endEMBL, V);
  595.     if (!V->allDone) {
  596.       while (!(EOFreader(V) ||
  597.          ((*V->s != '\0') && (strstr(V->s,"ID   ") == V->s))))
  598.       getline(V);
  599.     }
  600.     if (EOFreader(V)) V->allDone = true;
  601.   }
  602. }
  603.  
  604.  
  605.  
  606. Local boolean endZuker( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  607. {
  608.   *addend = false;
  609.   *ungetend= true;
  610.   return( *V->s == '(' );
  611. }
  612.  
  613. Local void readZuker(struct ReadSeqVars *V)
  614. {
  615.   /*! 1st string is Zuker's Fortran format */
  616.  
  617.   while (!V->allDone) {
  618.     getline(V);  /*s == "seqLen seqid string..."*/
  619.     strcpy(V->seqid, (V->s)+6);
  620.     readLoop(0, false, endZuker, V);
  621.     if (!V->allDone) {
  622.       while (!(EOFreader(V) |
  623.         ((*V->s != '\0') && (*V->s == '('))))
  624.           getline(V);
  625.       }
  626.     if (EOFreader(V)) V->allDone = true;
  627.   }
  628. }
  629.  
  630.  
  631.  
  632. Local boolean endFitch( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  633. {
  634.   /* this is a somewhat shaky end,
  635.     1st char of line is non-blank for seqq. title
  636.   */
  637.   *addend = false;
  638.   *ungetend= true;
  639.   return( *V->s != ' ' );
  640. }
  641.  
  642. Local void readFitch(struct ReadSeqVars *V)
  643. {
  644.   boolean first;
  645.  
  646.   first = true;
  647.   while (!V->allDone) {
  648.     if (!first) strcpy(V->seqid, V->s);
  649.     readLoop(0, first, endFitch, V);
  650.     if (EOFreader(V)) V->allDone = true;
  651.     first = false;
  652.     }
  653. }
  654.  
  655.  
  656. Local void readPlain(struct ReadSeqVars *V)
  657. {
  658.   V->nseq++;
  659.   V->addit = (V->choice > 0);
  660.   if (V->addit) V->seqlen = 0;
  661.   addseq(V->seqid, V);   /*from above..*/
  662.   if (V->fname!=NULL) sprintf(V->seqid, "%s  [Unknown form]", V->fname);
  663.   else sprintf(V->seqid, "  [Unknown form]");
  664.   do {
  665.     addseq(V->s, V);
  666.     V->done = EOFreader(V);
  667.     getline(V);
  668.   } while (!V->done);
  669.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  670.   V->allDone = true;
  671. }
  672.  
  673.  
  674. Local void readUWGCG(struct ReadSeqVars *V)
  675. {
  676. /*
  677. 10nov91: Reading GCG files casued duplication of last line when
  678.          EOF followed that line !!!
  679.     fix: getline now sets *V->s = 0
  680. */
  681.   char  *si;
  682.  
  683.   V->nseq++;
  684.   V->addit = (V->choice > 0);
  685.   if (V->addit) V->seqlen = 0;
  686.   strcpy(V->seqid, V->s);
  687.   /*writeseq: "    %s  Length: %d  (today)  Check: %d  ..\n" */
  688.   /*drop above or ".." from id*/
  689.   if ((si = strstr(V->seqid,"  Length: ")) != 0) *si = 0;
  690.   else if ((si = strstr(V->seqid,"..")) != 0) *si = 0;
  691.   do {
  692.     V->done = EOFreader(V);
  693.     getline(V);
  694.     if (!V->done) addseq((V->s), V);
  695.   } while (!V->done);
  696.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  697.   V->allDone = true;
  698. }
  699.  
  700.  
  701. Local void readOlsen(struct ReadSeqVars *V)
  702. { /* G. Olsen /print output from multiple sequence editor */
  703.  
  704.   char    *si, *sj, *sk, *sm, sid[40], snum[20];
  705.   boolean indata = false;
  706.   short snumlen = 0;
  707.  
  708.   V->addit = (V->choice > 0);
  709.   if (V->addit) V->seqlen = 0;
  710.   REWINDreader(V); 
  711.     V->nseq= 0;
  712.   do {
  713.     getline(V);
  714.     V->done = EOFreader(V);
  715.  
  716.     if (V->done && !(*V->s)) break;
  717.     else if (indata) {
  718.       if ( (si= strstr(V->s, sid))
  719.         && ( (sm= strstr(V->s, snum)) != NULL ) && (sm < si - snumlen) ) {
  720.         /* && (strstr(V->s, snum) == si - snumlen - 1) ) */
  721.  
  722.         /* Spaces are valid alignment data !! */
  723. /* 17Oct91: Error, the left margin is 21 not 22! */
  724. /* dropped some nucs up to now -- my example file was right shifted ! */
  725. /* variable right id margin, drop id-2 spaces at end */
  726. /*
  727.   VMS CC COMPILER (VAXC031) mess up:
  728.   -- Index of 21 is chopping 1st nuc on VMS systems Only!
  729.   Byte-for-byte same ame rnasep.olsen sequence file !
  730. */
  731.  
  732.         /* si = (V->s)+21; < was this before VMS CC wasted my time */
  733.         si += 10;  /* use strstr index plus offset to outfox VMS CC bug */
  734.  
  735.         if ((sk = strstr(si, sid))!=0) *(sk-2) = 0;
  736.         for (sk = si; *sk != 0; sk++) {
  737.            if (*sk == ' ') *sk = '.';
  738.            /* 18aug92: !! some olsen masks are NUMBERS !! which addseq eats */
  739.            else if (isdigit(*sk)) *sk= nonummask[*sk - '0'];
  740.            }
  741.  
  742.         addseq(si, V);
  743.         }
  744.       }
  745.  
  746.     else if ((sk = strstr(V->s, "): "))!=0) {  /* seqq info header line */
  747.   /* 18aug92: correct for diff seqs w/ same name -- use number, e.g. */
  748.   /*   3 (Agr.tume):  agrobacterium.prna  18-JUN-1987 16:12 */
  749.   /* 328 (Agr.tume):  agrobacterium.prna XYZ  19-DEC-1992   */
  750.       (V->nseq)++;
  751.       si = 1 + strchr(V->s,'(');
  752.       *sk = ' ';
  753.       if (V->choice == kListSequences) addinfo( si, V);
  754.       else if (V->nseq == V->choice) {
  755.         strcpy(V->seqid, si);
  756.         sj = strchr(V->seqid, ':');
  757.         while (*(--sj) == ' ') ;
  758.         while (--sj != V->seqid) { if (*sj == ' ') *sj = '_'; }
  759.  
  760.         *sk = 0;
  761.         while (*(--sk) == ' ') *sk = 0;
  762.         strcpy(sid, si);
  763.  
  764.         si= V->s;
  765.         while ((*si <= ' ') && (*si != 0)) si++;
  766.         snumlen=0;
  767.         while (si[snumlen] > ' ' && snumlen<20)
  768.          { snum[snumlen]= si[snumlen]; snumlen++; }
  769.         snum[snumlen]= 0;
  770.         }
  771.  
  772.       }
  773.  
  774.     else if (strstr(V->s,"identity:   Data:")) {
  775.       indata = true;
  776.       if (V->choice == kListSequences) V->done = true;
  777.       }
  778.  
  779.   } while (!V->done);
  780.  
  781.   V->allDone = true;
  782. } /*readOlsen*/
  783.  
  784.  
  785. Local void readMSF(struct ReadSeqVars *V)
  786. { /* gcg's MSF, mult. sequence format, interleaved ! */
  787.  
  788.   char    *si, *sj, sid[128];
  789.   boolean indata = false;
  790.   int     atseq= 0, iline= 0;
  791.  
  792.   V->addit = (V->choice > 0);
  793.   if (V->addit) V->seqlen = 0;
  794.   REWINDreader(V); 
  795.     V->nseq= 0;
  796.   do {
  797.     getline(V);
  798.     V->done = EOFreader(V);
  799.  
  800.     if (V->done && !(*V->s)) break;
  801.     else if (indata) {
  802.       /*somename  ...gpvedai .......t.. aaigr..vad tvgtgptnse aipaltaaet */
  803.       /*       E  gvenae.kgv tentna.tad fvaqpvylpe .nqt...... kv.affynrs */
  804.  
  805.       si= V->s;
  806.       skipwhitespace(si);
  807.       /* for (sj= si; isalnum(*sj); sj++) ; bug -- delwiche uses "-", "_" and others in names*/
  808.       for (sj= si; *sj > ' '; sj++) ;
  809.       *sj= 0;
  810.       if ( *si ) {
  811.         if ( (0==strcmp(si, sid)) ) {
  812.           addseq(sj+1, V);
  813.           }
  814.         iline++;
  815.         }
  816.       }
  817.  
  818.     else if (NULL != (si = strstr(V->s, "Name: "))) {  /* seqq info header line */
  819.       /* Name: somename      Len:   100  Check: 7009  Weight:  1.00 */
  820.  
  821.       (V->nseq)++;
  822.       si += 6;
  823.       if (V->choice == kListSequences) addinfo( si, V);
  824.       else if (V->nseq == V->choice) {
  825.         strcpy(V->seqid, si);
  826.         si = V->seqid;
  827.         skipwhitespace(si);
  828.         /* for (sj= si; isalnum(*sj); sj++) ; -- bug */
  829.         for (sj= si; *sj > ' '; sj++) ;
  830.         *sj= 0;
  831.         strcpy(sid, si);
  832.         }
  833.       }
  834.  
  835.     else if ( strstr(V->s,"//") /*== V->s*/ )  {
  836.       indata = true;
  837.       iline= 0;
  838.       if (V->choice == kListSequences) V->done = true;
  839.       }
  840.  
  841.   } while (!V->done);
  842.  
  843.  
  844.   V->allDone = true;
  845. } /*readMSF*/
  846.  
  847.  
  848.  
  849. Local void readPAUPinterleaved(struct ReadSeqVars *V)
  850. { /* PAUP mult. sequence format, interleaved or sequential! */
  851.  
  852.   char    *si, *sj, *send, sid[40], sid1[40], saveseq[255];
  853.   boolean first = true, indata = false, domatch;
  854.   int     atseq= 0, iline= 0, ifmc, saveseqlen=0;
  855.  
  856. #define fixmatchchar(s) { \
  857.   for (ifmc=0; ifmc<saveseqlen; ifmc++) \
  858.     if (s[ifmc] == V->matchchar) s[ifmc]= saveseq[ifmc]; }
  859.  
  860.   V->addit = (V->choice > 0);
  861.   V->seqlencount = 0;
  862.   if (V->addit) V->seqlen = 0;
  863.   /* rewind(V->f); V->nseq= 0;  << do in caller !*/
  864.   indata= true; /* call here after we find "matrix" */
  865.   domatch= (V->matchchar > 0);
  866.  
  867.   do {
  868.     getline(V);
  869.     V->done = EOFreader(V);
  870.  
  871.     if (V->done && !(*V->s)) break;
  872.     else if (indata) {
  873.       /* [         1                    1                    1         ]*/
  874.       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  875.       /* chimp     ................a.t. .c.................a ..........*/
  876.       /* !! need to correct for V->matchchar */
  877.       si= V->s;
  878.       skipwhitespace(si);
  879.       if (strchr(si,';')) indata= false;
  880.  
  881.       if (isalnum(*si))  {
  882.         /* valid data line starts w/ a left-justified seqq name in columns [0..8] */
  883.         if (first) {
  884.           (V->nseq)++;
  885.           if (V->nseq >= V->topnseq) first= false;
  886.           for (sj = si; isalnum(*sj); sj++) ;
  887.           send= sj;
  888.           skipwhitespace(sj);
  889.           if (V->choice == kListSequences) {
  890.             *send= 0;
  891.             addinfo( si, V);
  892.             }
  893.           else if (V->nseq == V->choice) {
  894.             if (domatch) {
  895.               if (V->nseq == 1) { strcpy( saveseq, sj); saveseqlen= strlen(saveseq); }
  896.               else fixmatchchar( sj);
  897.               }
  898.             addseq(sj, V);
  899.             *send= 0;
  900.             strcpy(V->seqid, si);
  901.             strcpy(sid, si);
  902.             if (V->nseq == 1) strcpy(sid1, sid);
  903.             }
  904.           }
  905.  
  906.         else if ( (strstr(si, sid) == si) ){
  907.           while (isalnum(*si)) si++;
  908.           skipwhitespace(si);
  909.           if (domatch) {
  910.             if (V->nseq == 1) { strcpy( saveseq, si); saveseqlen= strlen(saveseq); }
  911.             else fixmatchchar( si);
  912.             }
  913.           addseq(si, V);
  914.           }
  915.  
  916.         else if (domatch && (strstr(si, sid1) == si)) {
  917.           strcpy( saveseq, si);
  918.           saveseqlen= strlen(saveseq);
  919.           }
  920.  
  921.         iline++;
  922.         }
  923.       }
  924.  
  925.     else if ( strstr(V->s,"matrix") )  {
  926.       indata = true;
  927.       iline= 0;
  928.       if (V->choice == kListSequences) V->done = true;
  929.       }
  930.  
  931.   } while (!V->done);
  932.  
  933.   V->allDone = true;
  934. } /*readPAUPinterleaved*/
  935.  
  936.  
  937.  
  938. Local void readPAUPsequential(struct ReadSeqVars *V)
  939. { /* PAUP mult. sequence format, interleaved or sequential! */
  940.   char    *si, *sj;
  941.   boolean atname = true, indata = false;
  942.  
  943.   V->addit = (V->choice > 0);
  944.   if (V->addit) V->seqlen = 0;
  945.   V->seqlencount = 0;
  946.   /* REWINDreader(V); V->nseq= 0;  << do in caller !*/
  947.   indata= true; /* call here after we find "matrix" */
  948.   do {
  949.     getline(V);
  950.     V->done = EOFreader(V);
  951.  
  952.     if (V->done && !(*V->s)) break;
  953.     else if (indata) {
  954.       /* [         1                    1                    1         ]*/
  955.       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  956.       /*           aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  957.       /* chimp     ................a.t. .c.................a ..........*/
  958.       /*           ................a.t. .c.................a ..........*/
  959.  
  960.       si= V->s;
  961.       skipwhitespace(si);
  962.       if (strchr(si,';')) indata= false;
  963.       if (isalnum(*si))  {
  964.         /* valid data line starts w/ a left-justified seqq name in columns [0..8] */
  965.         if (atname) {
  966.           (V->nseq)++;
  967.           V->seqlencount = 0;
  968.           atname= false;
  969.           sj= si+1;
  970.           while (isalnum(*sj)) sj++;
  971.           if (V->choice == kListSequences) {
  972.             /* !! we must count bases to know when topseqlen is reached ! */
  973.             countseq(sj, V);
  974.             if (V->seqlencount >= V->topseqlen) atname= true;
  975.             *sj= 0;
  976.             addinfo( si, V);
  977.             }
  978.           else if (V->nseq == V->choice) {
  979.             addseq(sj, V);
  980.             V->seqlencount= V->seqlen;
  981.             if (V->seqlencount >= V->topseqlen) atname= true;
  982.             *sj= 0;
  983.             strcpy(V->seqid, si);
  984.             }
  985.           else {
  986.             countseq(sj, V);
  987.             if (V->seqlencount >= V->topseqlen) atname= true;
  988.             }
  989.           }
  990.  
  991.         else if (V->nseq == V->choice) {
  992.           addseq(V->s, V);
  993.           V->seqlencount= V->seqlen;
  994.           if (V->seqlencount >= V->topseqlen) atname= true;
  995.           }
  996.         else {
  997.           countseq(V->s, V);
  998.           if (V->seqlencount >= V->topseqlen) atname= true;
  999.           }
  1000.         }
  1001.       }
  1002.  
  1003.     else if ( strstr(V->s,"matrix") )  {
  1004.       indata = true;
  1005.       atname= true;
  1006.       if (V->choice == kListSequences) V->done = true;
  1007.       }
  1008.  
  1009.   } while (!V->done);
  1010.  
  1011.   V->allDone = true;
  1012. } /*readPAUPsequential*/
  1013.  
  1014.  
  1015. Local void readPhylipInterleaved(struct ReadSeqVars *V)
  1016. {
  1017.   char    *si, *sj;
  1018.   boolean first = true;
  1019.   int     iline= 0;
  1020.  
  1021.   V->addit = (V->choice > 0);
  1022.   if (V->addit) V->seqlen = 0;
  1023.   V->seqlencount = 0;
  1024.   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); << topnseq == 0 !!! bad scan !! */
  1025.   si= V->s;
  1026.   skipwhitespace(si);
  1027.   V->topnseq= atoi(si);
  1028.   while (isdigit(*si)) si++;
  1029.   skipwhitespace(si);
  1030.   V->topseqlen= atol(si);
  1031.   /* fprintf(stderr,"Phylip-ileaf: topnseq=%d  topseqlen=%d\n",V->topnseq, V->topseqlen); */
  1032.  
  1033.   do {
  1034.     getline(V);
  1035.     V->done = EOFreader(V);
  1036.  
  1037.     if (V->done && !(*V->s)) break;
  1038.     si= V->s;
  1039.     skipwhitespace(si);
  1040.     if (*si != 0) {
  1041.  
  1042.       if (first) {  /* collect seqq names + seqq, as fprintf(outf,"%-10s  ",seqname); */
  1043.         (V->nseq)++;
  1044.         if (V->nseq >= V->topnseq) first= false;
  1045.         sj= V->s+10;  /* past name, start of data */
  1046.         if (V->choice == kListSequences) {
  1047.           *sj= 0;
  1048.           addinfo( si, V);
  1049.           }
  1050.         else if (V->nseq == V->choice) {
  1051.           addseq(sj, V);
  1052.           *sj= 0;
  1053.           strcpy(V->seqid, si);
  1054.           }
  1055.         }
  1056.       else if ( iline % V->nseq == V->choice -1 ) {
  1057.         addseq(si, V);
  1058.         }
  1059.       iline++;
  1060.     }
  1061.   } while (!V->done);
  1062.  
  1063.   V->allDone = true;
  1064. } /*readPhylipInterleaved*/
  1065.  
  1066.  
  1067.  
  1068. Local boolean endPhylipSequential( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  1069. {
  1070.   boolean done;
  1071.     /* *addend = false; << need this true if EOF(file) on last dataline */
  1072.   *ungetend= false;
  1073.   countseq( V->s, V);
  1074.     done= (V->seqlencount >= V->topseqlen);
  1075.   *addend = !done;
  1076.     return done;
  1077. }
  1078.  
  1079. Local void readPhylipSequential(struct ReadSeqVars *V)
  1080. {
  1081.   short  i;
  1082.   char  *si;
  1083.   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); < ? bad sscan ? */
  1084.   si= V->s;
  1085.   skipwhitespace(si);
  1086.   V->topnseq= atoi(si);
  1087.   while (isdigit(*si)) si++;
  1088.   skipwhitespace(si);
  1089.   V->topseqlen= atol(si);
  1090.   getline(V);
  1091.   while (!V->allDone) {
  1092.     V->seqlencount= 0;
  1093.     strncpy(V->seqid, (V->s), 10);
  1094.     V->seqid[10]= 0;
  1095.     for (i=0; i<10 && V->s[i]; i++) V->s[i]= ' ';
  1096.     readLoop(0, true, endPhylipSequential, V);
  1097.     if (EOFreader(V)) V->allDone = true;
  1098.     }
  1099. }
  1100.  
  1101.  
  1102.  
  1103.  
  1104. Local void readSeqMain(
  1105.       struct ReadSeqVars *V,
  1106.       const long  skiplines,
  1107.       const short format)
  1108. {
  1109. #define tolowerstr(s) { long Itlwr, Ntlwr= strlen(s); \
  1110.   for (Itlwr=0; Itlwr<Ntlwr; Itlwr++) s[Itlwr]= to_lower(s[Itlwr]); }
  1111.  
  1112.   boolean gotuw;
  1113.   long l;
  1114.  
  1115.     /*fprintf(stderr,"readSeqMain &reader= %ld  &gFile= %ld\n", V->reader, gFile);*/
  1116.  
  1117.   V->linestart= 0;
  1118.   V->matchchar= 0;
  1119.   if (V->reader == NULL)
  1120.     V->err = eFileNotFound;
  1121.   else {
  1122.  
  1123.     for (l = skiplines; l > 0; l--) getline( V);
  1124.  
  1125.     do {
  1126.       getline( V);
  1127.       for (l= strlen(V->s); (l > 0) && (V->s[l] == ' '); l--) ;
  1128.     } while ((l == 0) && !EOFreader(V));
  1129.  
  1130.     if (EOFreader(V)) V->err = eNoData;
  1131.  
  1132.     else switch (format) {
  1133.       case kPlain : readPlain(V); break;
  1134.       case kIG    : readIG(V); break;
  1135.       case kStrider: readStrider(V); break;
  1136.       case kGenBank: readGenBank(V); break;
  1137.       case kPIR   : readPIR(V); break;
  1138.       case kNBRF  : readNBRF(V); break;
  1139.       case kPearson: readPearson(V); break;
  1140.       case kEMBL  : readEMBL(V); break;
  1141.       case kZuker : readZuker(V); break;
  1142.       case kOlsen : readOlsen(V); break;
  1143.       case kMSF   : 
  1144.           V->isseqchar = isGCGSeqChar;
  1145.                 readMSF(V); 
  1146.                 break;
  1147.  
  1148.       case kPAUP    : {
  1149.         boolean done= false;
  1150.         boolean interleaved= false;
  1151.         char  *cp;
  1152.         /* REWINDreader(V); V->nseq= 0; ?? assume it is at top ?? skiplines ... */
  1153.         while (!done) {
  1154.           getline( V);
  1155.           tolowerstr( V->s);
  1156.           if (strstr( V->s, "matrix")) done= true;
  1157.           if (strstr( V->s, "interleav")) interleaved= true;
  1158.           if (NULL != (cp=strstr( V->s, "ntax=")) )  V->topnseq= atoi(cp+5);
  1159.           if (NULL != (cp=strstr( V->s, "nchar=")) )  V->topseqlen= atoi(cp+6);
  1160.           if (NULL != (cp=strstr( V->s, "matchchar=")) )  {
  1161.             cp += 10;
  1162.             if (*cp=='\'') cp++;
  1163.             else if (*cp=='"') cp++;
  1164.             V->matchchar= *cp;
  1165.             }
  1166.           }
  1167.         if (interleaved) readPAUPinterleaved(V);
  1168.         else readPAUPsequential(V);
  1169.         }
  1170.         break;
  1171.  
  1172.       /* kPhylip: ! can't determine in middle of file which type it is...*/
  1173.       /* test for interleave or sequential and use Phylip4(ileave) or Phylip2 */
  1174.       case kPhylip2:
  1175.         readPhylipSequential(V);
  1176.         break;
  1177.       case kPhylip4: /* == kPhylip3 */
  1178.         readPhylipInterleaved(V);
  1179.         break;
  1180.  
  1181.       default:
  1182.         V->err = eUnknownFormat;
  1183.         break;
  1184.  
  1185.       case kFitch :
  1186.         strcpy(V->seqid, V->s); getline(V);
  1187.         readFitch(V);
  1188.         break;
  1189.  
  1190.       case kGCG:
  1191.         V->isseqchar = isGCGSeqChar;
  1192.         do {
  1193.           gotuw = (strstr(V->s,"..") != NULL);
  1194.           if (gotuw) readUWGCG(V);
  1195.           getline(V);
  1196.         } while (!(EOFreader(V) || V->allDone));
  1197.         break;
  1198.       }
  1199.     }
  1200.  
  1201.   V->filestart= false;
  1202.   V->seq[V->seqlen] = 0; /* stick a string terminator on it */
  1203.     REALLOC(V->seq, V->seqlen+1); /* and resize it down to true length */
  1204. }
  1205.  
  1206.  
  1207. char *readSeqFp(
  1208.       const short whichEntry,  /* index to sequence in file */
  1209.       FILE  *fp,   /* pointer to open seqq file */
  1210.       const long  skiplines,
  1211.       const short format,      /* sequence file format */
  1212.       long  *seqlen,     /* return seqq size */
  1213.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1214.       short *error,      /* return error */
  1215.       char  *seqid)      /* return seqq name/info */
  1216. {
  1217.   struct ReadSeqVars V;
  1218.  
  1219.     gFile= fp; gLinestart = 0;
  1220.     V.reader = readWriteFromFile;
  1221.     
  1222.   if (format < kMinFormat || format > kMaxFormat) {
  1223.     *error = eUnknownFormat;
  1224.     *seqlen = 0;
  1225.     return NULL;
  1226.     }
  1227.  
  1228.   V.choice = whichEntry;
  1229.   V.fname  = NULL;  /* don't know */
  1230.   V.seq    = (char*) MALLOC(kStartLength+1);
  1231.   V.maxseq = kStartLength;
  1232.   V.memstep= 0;
  1233.   aStartLength= kStartLength;
  1234.   V.seqlen = 0;
  1235.   V.seqid  = seqid;
  1236.  
  1237.   V.f = fp;
  1238.   V.filestart= (ftell( fp) == 0); 
  1239.   /* !! in sequential read, must remove current seqq position from choice/whichEntry counter !! ... */
  1240.   if (V.filestart)  V.nseq = 0;
  1241.   else V.nseq= *nseq;  /* track where we are in file...*/
  1242.  
  1243.   *V.seqid = '\0';
  1244.   V.err = 0;
  1245.   V.nseq = 0;
  1246.   V.isseqchar = isSeqChar;
  1247.   if (V.choice == kListSequences) ; /* leave as is */
  1248.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1249.   V.addit = (V.choice > 0);
  1250.   V.allDone = false;
  1251.  
  1252. #if 0
  1253.      gLinestart= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1254.   (*V.reader)((char*)NULL, 0L, kRSFile_End);
  1255.      V.maxseq= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1256.   (*V.reader)((char*)NULL, gLinestart, kRSFile_Seek);
  1257. #endif
  1258.  
  1259.   readSeqMain(&V, skiplines, format);
  1260.  
  1261.   *error = V.err;
  1262.   *seqlen = V.seqlen;
  1263.   *nseq = V.nseq;
  1264.   return V.seq;
  1265. }
  1266.  
  1267. char *readSeq(
  1268.       const short whichEntry,  /* index to sequence in file */
  1269.       const char  *filename,   /* file name */
  1270.       const long  skiplines,
  1271.       const short format,      /* sequence file format */
  1272.       long  *seqlen,     /* return seqq size */
  1273.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1274.       short *error,      /* return error */
  1275.       char  *seqid)      /* return seqq name/info */
  1276. {
  1277.   struct ReadSeqVars V;
  1278.  
  1279.   if (format < kMinFormat || format > kMaxFormat) {
  1280.     *error = eUnknownFormat;
  1281.     *seqlen = 0;
  1282.     return NULL;
  1283.     }
  1284.  
  1285.   V.choice = whichEntry;
  1286.   V.fname  = filename;  /* don't need to copy string, just ptr to it */
  1287.   V.seq    = (char*) MALLOC( kStartLength+1);
  1288.   V.maxseq = kStartLength;
  1289.   V.memstep= 0;
  1290.   aStartLength= kStartLength;
  1291.   V.seqlen = 0;
  1292.   V.seqid  = seqid;
  1293.  
  1294.   V.f = NULL;
  1295.   *V.seqid = '\0';
  1296.   V.err = 0;
  1297.   V.nseq = 0;
  1298.   V.isseqchar = isSeqChar;
  1299.   if (V.choice == kListSequences) ; /* leave as is */
  1300.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1301.   V.addit = (V.choice > 0);
  1302.   V.allDone = false;
  1303.  
  1304.   V.f = fopen(V.fname, "r");
  1305.   V.filestart= true;
  1306.     gFile= V.f; gLinestart = 0;
  1307.     V.reader = readWriteFromFile;
  1308. #if 0
  1309.      gLinestart= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1310.   (*V.reader)((char*)NULL, 0L, kRSFile_End);
  1311.      V.maxseq= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1312.   (*V.reader)((char*)NULL, gLinestart, kRSFile_Seek);
  1313. #endif
  1314.  
  1315.   readSeqMain(&V, skiplines, format);
  1316.  
  1317.   if (V.f != NULL) fclose(V.f);
  1318.   *error = V.err;
  1319.   *seqlen = V.seqlen;
  1320.   *nseq = V.nseq;
  1321.   return V.seq;
  1322. }
  1323.  
  1324.  
  1325.  
  1326. char *readSeqCall(
  1327.       const short whichEntry,  /* index to sequence in file */
  1328.         ReadWriteProc reader,
  1329.       const long  skiplines,
  1330.       const short format,      /* sequence file format */
  1331.       long  *seqlen,     /* return seqq size */
  1332.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1333.       short *error,      /* return error */
  1334.       char  *seqid)      /* return seqq name/info */
  1335. {
  1336.   struct ReadSeqVars V;
  1337.  
  1338.   if (format < kMinFormat || format > kMaxFormat) {
  1339.     *error = eUnknownFormat;
  1340.     *seqlen = 0;
  1341.     return NULL;
  1342.     }
  1343.  
  1344.   V.choice = whichEntry;
  1345.   V.fname  = NULL;  /* don't know */
  1346.   V.seq    = (char*) MALLOC( kStartLength+1);
  1347.   V.maxseq = kStartLength;
  1348.   V.memstep= 0;
  1349.   aStartLength= kStartLength;
  1350.   V.seqlen = 0;
  1351.   V.seqid  = seqid;
  1352.   *V.seqid = '\0';
  1353.  
  1354.   V.f = NULL;
  1355.   V.err = 0;
  1356.   V.nseq = 0;
  1357.   V.isseqchar = isSeqChar;
  1358.   if (V.choice == kListSequences) ; /* leave as is */
  1359.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1360.   V.addit = (V.choice > 0);
  1361.   V.allDone = false;
  1362.  
  1363.     gFile= NULL; gLinestart = 0;
  1364.     V.reader = reader;
  1365.      V.filestart= (0 == (*reader)((char*)NULL, 0L, kRSFile_Tell));
  1366. #if 0
  1367.      gLinestart= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1368.   (*V.reader)((char*)NULL, 0L, kRSFile_End);
  1369.      V.maxseq= (*V.reader)((char*)NULL, 0L, kRSFile_Tell);
  1370.   (*V.reader)((char*)NULL, gLinestart, kRSFile_Seek);
  1371. #endif
  1372.  
  1373.   readSeqMain(&V, skiplines, format);
  1374.  
  1375.   *error = V.err;
  1376.   *seqlen = V.seqlen;
  1377.   *nseq = V.nseq;
  1378.   return V.seq;
  1379. }
  1380.  
  1381.  
  1382.  
  1383. char *listSeqs(
  1384.       const char  *filename,   /* file name */
  1385.       const long skiplines,
  1386.       const short format,      /* sequence file format */
  1387.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1388.       short *error)      /* return error */
  1389. {
  1390.   char  seqid[256];
  1391.   long  seqlen;
  1392.  
  1393.     /*fprintf(stderr,"listSeqs filename= %s  format= %d  nseq= %d\n", filename, format, *nseq);*/
  1394.   return readSeq( kListSequences, filename, skiplines, format,
  1395.                   &seqlen, nseq, error, seqid);
  1396. }
  1397.  
  1398. char *listSeqsCall(
  1399.       ReadWriteProc    reader,
  1400.       const long skiplines,
  1401.       const short format,      /* sequence file format */
  1402.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1403.       short *error)      /* return error */
  1404. {
  1405.   char  seqid[256];
  1406.   long  seqlen;
  1407.  
  1408.   return readSeqCall( kListSequences, reader, skiplines, format,
  1409.                   &seqlen, nseq, error, seqid);
  1410. }
  1411.  
  1412.  
  1413.  
  1414.  
  1415. short seqFileFormat(    /* return sequence format number, see ureadseq.h */
  1416.     const char *filename,
  1417.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1418.     short *error)       /* return any error value or 0 */
  1419. {
  1420.   FILE      *fseq;
  1421.   short      format;
  1422.  
  1423.   fseq  = fopen(filename, "r");
  1424.   format= seqFileFormatFp( fseq, skiplines, error);
  1425.   if (fseq!=NULL) fclose(fseq);
  1426.   return format;
  1427. }
  1428.  
  1429. short seqFileFormatFp(
  1430.     FILE *fseq,
  1431.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1432.     short *error)       /* return any error value or 0 */
  1433. {
  1434.     gFile = fseq; gLinestart = 0;
  1435.     return seqFileFormatCall( readWriteFromFile, skiplines, error);
  1436. }
  1437.  
  1438.         
  1439.  
  1440. short seqFileFormatCall(
  1441.     ReadWriteProc reader,
  1442.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1443.     short *error)       /* return any error value or 0 */
  1444. {
  1445.   boolean   foundDNA= false, foundIG= false, foundStrider= false,
  1446.             foundGB= false, foundPIR= false, foundEMBL= false, foundNBRF= false,
  1447.             foundPearson= false, foundFitch= false, foundPhylip= false, foundZuker= false,
  1448.             gotolsen= false, gotpaup = false, gotasn1 = false, gotuw= false, gotMSF= false,
  1449.             isfitch= false,  isphylip= false, done= false;
  1450.   short     format= kUnknown;
  1451.   int       nlines= 0, k, splen= 0, otherlines= 0, aminolines= 0, dnalines= 0;
  1452.   char      sp[256];
  1453.   long      linestart=0;
  1454.   int         maxlines2check=500;
  1455.  
  1456. #define ReadOneLine(sp)   \
  1457.          { done |= (boolean)((*reader)((char*)NULL, 0L, kRSFile_Eof)); \
  1458.          callreadline( reader, sp);  \
  1459.     if (!done) { splen = strlen(sp); ++nlines; } }
  1460.  
  1461. /*  { done |= feof(fp); \ */
  1462. /*      readline( fp, sp, &linestart); \ */
  1463.  
  1464.     /*fprintf(stderr,"seqFileFormatCall &reader= %ld  &gFile= %ld\n", reader, gFile); */
  1465.  
  1466.   *skiplines = 0;
  1467.   *error = 0;
  1468.   if (reader == NULL) { *error = eFileNotFound;  return kNoformat; }
  1469.  
  1470.   while ( !done ) {
  1471.     ReadOneLine(sp);
  1472.  
  1473.     /* check for mailer head & skip past if found */
  1474.     if (nlines < 4 && !done) {
  1475.       if ((strstr(sp,"From ") == sp) || (strstr(sp,"Received:") == sp)) {
  1476.         do {
  1477.           /* skip all lines until find one blank line */
  1478.           ReadOneLine(sp);
  1479.           if (!done) for (k=0; (k<splen) && (sp[k]==' '); k++) ;
  1480.           } while ((!done) && (k < splen));
  1481.         *skiplines = nlines; /* !? do we want #lines or #bytes ?? */
  1482.         }
  1483.       }
  1484.  
  1485.     if (sp==NULL || *sp==0)
  1486.       ; /* nada */
  1487.  
  1488.     /* high probability identities: */
  1489.  
  1490.     else if ( strstr(sp,"MSF:") && strstr(sp,"Type:") && strstr(sp,"Check:") )
  1491.       gotMSF= true;
  1492.  
  1493.     else if ((strstr(sp,"..") != NULL) && (strstr(sp,"Check:") != NULL))
  1494.       gotuw= true;
  1495.  
  1496.     else if (strstr(sp,"identity:   Data:") != NULL)
  1497.       gotolsen= true;
  1498.  
  1499.     else if ( strstr(sp,"::=") &&
  1500.       (strstr(sp,"Bioseq") ||       /* Bioseq or Bioseq-set */
  1501.        strstr(sp,"Seq-entry") ||
  1502.        strstr(sp,"Seq-submit") ) )  /* can we read submit format? */
  1503.           gotasn1= true;
  1504.  
  1505.     else if ( strstr(sp,"#NEXUS") == sp )
  1506.       gotpaup= true;
  1507.  
  1508.     /* uncertain identities: */
  1509.  
  1510.     else if (*sp ==';') {
  1511.       if (strstr(sp,"Strider") !=NULL) foundStrider= true;
  1512.       else foundIG= true;
  1513.       }
  1514.  
  1515.     else if (strstr(sp,"LOCUS") == sp)
  1516.       foundGB= true;
  1517.     else if (strstr(sp,"ORIGIN") == sp)
  1518.       foundGB= true;
  1519.  
  1520.     else if (strstr(sp,"ENTRY   ") == sp)  /* ? also (strcmp(sp,"\\\\\\")==0) */
  1521.       foundPIR= true;
  1522.     else if (strstr(sp,"SEQUENCE") == sp)
  1523.       foundPIR= true;
  1524.  
  1525.     else if (*sp == '>') {
  1526.       if (sp[3] == ';') foundNBRF= true;
  1527.       else foundPearson= true;
  1528.       }
  1529.  
  1530.     else if (strstr(sp,"ID   ") == sp)
  1531.       foundEMBL= true;
  1532.     else if (strstr(sp,"SQ   ") == sp)
  1533.       foundEMBL= true;
  1534.  
  1535.     else if (*sp == '(')
  1536.       foundZuker= true;
  1537.  
  1538.     else {
  1539.       if (nlines - *skiplines == 1) {
  1540.         int ispp= 0, ilen= 0;
  1541.         sscanf( sp, "%d%d", &ispp, &ilen);
  1542.         if (ispp > 0 && ilen > 0) isphylip= true;
  1543.         }
  1544.       else if (isphylip && nlines - *skiplines == 2) {
  1545.         int  tseq;
  1546.         tseq= getseqtype(sp+10, strlen(sp+10));
  1547.         if ( isalpha(*sp)     /* 1st letter in 2nd line must be of a name */
  1548.          && (tseq != kOtherSeq))  /* sequence section must be okay */
  1549.             foundPhylip= true;
  1550.         }
  1551.  
  1552.       for (k=0, isfitch= true; isfitch && (k < splen); k++) {
  1553.         if (k % 4 == 0) isfitch &= (sp[k] == ' ');
  1554.         else isfitch &= (sp[k] != ' ');
  1555.         }
  1556.       if (isfitch && (splen > 20)) foundFitch= true;
  1557.  
  1558.       /* kRNA && kDNA are fairly certain...*/
  1559.       switch (getseqtype( sp, splen)) {
  1560.         case kOtherSeq: otherlines++; break;
  1561.         case kAmino   : if (splen>20) aminolines++; break;
  1562.         case kDNA     :
  1563.         case kRNA     : if (splen>20) dnalines++; break;
  1564.         case kNucleic : break; /* not much info ? */
  1565.         }
  1566.  
  1567.       }
  1568.  
  1569.                     /* pretty certain */
  1570.     if (gotolsen) {
  1571.       format= kOlsen;
  1572.       done= true;
  1573.       }
  1574.     else if (gotMSF) {
  1575.       format= kMSF;
  1576.       done= true;
  1577.       }
  1578.     else if (gotasn1) {
  1579.       /* !! we need to look further and return  kASNseqentry | kASNseqset */
  1580.       /*
  1581.         seqentry key is Seq-entry ::=
  1582.         seqset key is Bioseq-set ::=
  1583.         ?? can't read these yet w/ ncbi tools ??
  1584.           Seq-submit ::=
  1585.           Bioseq ::=  << fails both bioseq-seqq and seq-entry parsers !
  1586.       */
  1587.       if (strstr(sp,"Bioseq-set")) format= kASNseqset;
  1588.       else if (strstr(sp,"Seq-entry")) format= kASNseqentry;
  1589.       else format= kASN1;  /* other form, we can't yet read... */
  1590.       done= true;
  1591.       }
  1592.     else if (gotpaup) {
  1593.       format= kPAUP;
  1594.       done= true;
  1595.       }
  1596.  
  1597.     else if (gotuw) {
  1598.       if (foundIG) format= kIG;  /* a TOIG file from GCG for certain */
  1599.       else format= kGCG;
  1600.       done= true;
  1601.       }
  1602.  
  1603.     else if ((dnalines > 1) || done || (nlines > maxlines2check)) {
  1604.           /* decide on most likely format */
  1605.           /* multichar idents: */
  1606.       if (foundStrider) format= kStrider;
  1607.       else if (foundGB) format= kGenBank;
  1608.       else if (foundPIR) format= kPIR;
  1609.       else if (foundEMBL) format= kEMBL;
  1610.       else if (foundNBRF) format= kNBRF;
  1611.           /* single char idents: */
  1612.       else if (foundIG) format= kIG;
  1613.       else if (foundPearson) format= kPearson;
  1614.       else if (foundZuker) format= kZuker;
  1615.           /* digit ident: */
  1616.       else if (foundPhylip) format= kPhylip;
  1617.           /* spacing ident: */
  1618.       else if (foundFitch) format= kFitch;
  1619.           /* no format chars: */
  1620.       else if (otherlines > 0) format= kUnknown;
  1621.       else if (dnalines > 1) format= kPlain;
  1622.       else if (aminolines > 1) format= kPlain;
  1623.       else format= kUnknown;
  1624.  
  1625.       done= true;
  1626.       }
  1627.  
  1628.     /* need this for possible long header in olsen format */
  1629.      else if (strstr(sp,"): ") != NULL)
  1630.        maxlines2check++;
  1631.     }
  1632.  
  1633.   if (format == kPhylip) {
  1634.     /* check for interleaved or sequential -- really messy */
  1635.     int tname, tseq;
  1636.     long i, j, nspp= 0, nlen= 0, ilen, leaf= 0, seq= 0;
  1637.     char  *ps;
  1638.  
  1639.         (void) (*reader)((char*)NULL, 0L, kRSFile_Rewind);
  1640.         /*rewind(fp);*/
  1641.     for (i=0; i < *skiplines; i++) ReadOneLine(sp);
  1642.     nlines= 0;
  1643.     ReadOneLine(sp);
  1644.     sscanf( sp, "%d%d", &nspp, &nlen);
  1645.     ReadOneLine(sp); /* 1st seqq line */
  1646.     for (ps= sp+10, ilen=0; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1647.  
  1648.     for (i= 1; i<nspp; i++) {
  1649.       ReadOneLine(sp);
  1650.  
  1651.       tseq= getseqtype(sp+10, strlen(sp+10));
  1652.       tname= getseqtype(sp, 10);
  1653.       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++)  ;
  1654.       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1655.  
  1656.       /* find probable interleaf or sequential ... */
  1657.       if (j>=9) seq += 10; /* pretty certain not ileaf */
  1658.       else {
  1659.         if (tseq != tname) leaf++; else seq++;
  1660.         if (tname == kDNA || tname == kRNA) seq++; else leaf++;
  1661.         }
  1662.  
  1663.       if (ilen <= nlen && j<9) {
  1664.         if (tname == kOtherSeq) leaf += 10;
  1665.         else if (tname == kAmino || tname == kDNA || tname == kRNA) seq++; else leaf++;
  1666.         }
  1667.       else if (ilen > nlen) {
  1668.         ilen= 0;
  1669.         }
  1670.       }
  1671.     for ( nspp *= 2 ; i<nspp; i++) {  /* this should be only bases if interleaf */
  1672.       ReadOneLine(sp);
  1673.  
  1674.       tseq= getseqtype(sp+10, strlen(sp+10));
  1675.       tname= getseqtype(sp, 10);
  1676.       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1677.       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++)  ;
  1678.       if (j<9) {
  1679.         if (tname == kOtherSeq) seq += 10;
  1680.         if (tseq != tname) seq++; else leaf++;
  1681.         if (tname == kDNA || tname == kRNA) leaf++; else seq++;
  1682.         }
  1683.       if (ilen > nlen) {
  1684.         if (j>9) leaf += 10; /* must be a name here for sequent */
  1685.         else if (tname == kOtherSeq) seq += 10;
  1686.         ilen= 0;
  1687.         }
  1688.       }
  1689.  
  1690.     if (leaf > seq) format= kPhylip4;
  1691.     else format= kPhylip2;
  1692.     }
  1693.  
  1694.     /*fprintf(stderr, "seqFileFormatCall &reader= %ld  &gFile= %ld  format= %ld\n", reader, gFile, format);*/
  1695.   return(format);
  1696. #undef  ReadOneLine
  1697. } /* SeqFileFormat */
  1698.  
  1699.  
  1700.  
  1701.  
  1702. unsigned long GCGchecksum( const char *seq, const long seqlen, unsigned long *checktotal)
  1703. /* GCGchecksum */
  1704. {
  1705.   register long  i, check = 0, count = 0;
  1706.     register char *sp = (char*) seq;
  1707.  
  1708.   for (i = 0; i < seqlen; i++) {
  1709.     count++;
  1710.         check += count * to_upper(*sp);
  1711.       sp++;
  1712.     if (count == 57) count = 0;
  1713.     }
  1714.   check %= 10000;
  1715.   *checktotal += check;
  1716.   *checktotal %= 10000;
  1717.   return check;
  1718. }
  1719.  
  1720. /* Table of CRC-32's of all single byte values (made by makecrc.c of ZIP source) */
  1721. const unsigned long crctab[] = {
  1722.   0x00000000L, 0x77073096L, 0xee0e612cL, 0x990951baL, 0x076dc419L,
  1723.   0x706af48fL, 0xe963a535L, 0x9e6495a3L, 0x0edb8832L, 0x79dcb8a4L,
  1724.   0xe0d5e91eL, 0x97d2d988L, 0x09b64c2bL, 0x7eb17cbdL, 0xe7b82d07L,
  1725.   0x90bf1d91L, 0x1db71064L, 0x6ab020f2L, 0xf3b97148L, 0x84be41deL,
  1726.   0x1adad47dL, 0x6ddde4ebL, 0xf4d4b551L, 0x83d385c7L, 0x136c9856L,
  1727.   0x646ba8c0L, 0xfd62f97aL, 0x8a65c9ecL, 0x14015c4fL, 0x63066cd9L,
  1728.   0xfa0f3d63L, 0x8d080df5L, 0x3b6e20c8L, 0x4c69105eL, 0xd56041e4L,
  1729.   0xa2677172L, 0x3c03e4d1L, 0x4b04d447L, 0xd20d85fdL, 0xa50ab56bL,
  1730.   0x35b5a8faL, 0x42b2986cL, 0xdbbbc9d6L, 0xacbcf940L, 0x32d86ce3L,
  1731.   0x45df5c75L, 0xdcd60dcfL, 0xabd13d59L, 0x26d930acL, 0x51de003aL,
  1732.   0xc8d75180L, 0xbfd06116L, 0x21b4f4b5L, 0x56b3c423L, 0xcfba9599L,
  1733.   0xb8bda50fL, 0x2802b89eL, 0x5f058808L, 0xc60cd9b2L, 0xb10be924L,
  1734.   0x2f6f7c87L, 0x58684c11L, 0xc1611dabL, 0xb6662d3dL, 0x76dc4190L,
  1735.   0x01db7106L, 0x98d220bcL, 0xefd5102aL, 0x71b18589L, 0x06b6b51fL,
  1736.   0x9fbfe4a5L, 0xe8b8d433L, 0x7807c9a2L, 0x0f00f934L, 0x9609a88eL,
  1737.   0xe10e9818L, 0x7f6a0dbbL, 0x086d3d2dL, 0x91646c97L, 0xe6635c01L,
  1738.   0x6b6b51f4L, 0x1c6c6162L, 0x856530d8L, 0xf262004eL, 0x6c0695edL,
  1739.   0x1b01a57bL, 0x8208f4c1L, 0xf50fc457L, 0x65b0d9c6L, 0x12b7e950L,
  1740.   0x8bbeb8eaL, 0xfcb9887cL, 0x62dd1ddfL, 0x15da2d49L, 0x8cd37cf3L,
  1741.   0xfbd44c65L, 0x4db26158L, 0x3ab551ceL, 0xa3bc0074L, 0xd4bb30e2L,
  1742.   0x4adfa541L, 0x3dd895d7L, 0xa4d1c46dL, 0xd3d6f4fbL, 0x4369e96aL,
  1743.   0x346ed9fcL, 0xad678846L, 0xda60b8d0L, 0x44042d73L, 0x33031de5L,
  1744.   0xaa0a4c5fL, 0xdd0d7cc9L, 0x5005713cL, 0x270241aaL, 0xbe0b1010L,
  1745.   0xc90c2086L, 0x5768b525L, 0x206f85b3L, 0xb966d409L, 0xce61e49fL,
  1746.   0x5edef90eL, 0x29d9c998L, 0xb0d09822L, 0xc7d7a8b4L, 0x59b33d17L,
  1747.   0x2eb40d81L, 0xb7bd5c3bL, 0xc0ba6cadL, 0xedb88320L, 0x9abfb3b6L,
  1748.   0x03b6e20cL, 0x74b1d29aL, 0xead54739L, 0x9dd277afL, 0x04db2615L,
  1749.   0x73dc1683L, 0xe3630b12L, 0x94643b84L, 0x0d6d6a3eL, 0x7a6a5aa8L,
  1750.   0xe40ecf0bL, 0x9309ff9dL, 0x0a00ae27L, 0x7d079eb1L, 0xf00f9344L,
  1751.   0x8708a3d2L, 0x1e01f268L, 0x6906c2feL, 0xf762575dL, 0x806567cbL,
  1752.   0x196c3671L, 0x6e6b06e7L, 0xfed41b76L, 0x89d32be0L, 0x10da7a5aL,
  1753.   0x67dd4accL, 0xf9b9df6fL, 0x8ebeeff9L, 0x17b7be43L, 0x60b08ed5L,
  1754.   0xd6d6a3e8L, 0xa1d1937eL, 0x38d8c2c4L, 0x4fdff252L, 0xd1bb67f1L,
  1755.   0xa6bc5767L, 0x3fb506ddL, 0x48b2364bL, 0xd80d2bdaL, 0xaf0a1b4cL,
  1756.   0x36034af6L, 0x41047a60L, 0xdf60efc3L, 0xa867df55L, 0x316e8eefL,
  1757.   0x4669be79L, 0xcb61b38cL, 0xbc66831aL, 0x256fd2a0L, 0x5268e236L,
  1758.   0xcc0c7795L, 0xbb0b4703L, 0x220216b9L, 0x5505262fL, 0xc5ba3bbeL,
  1759.   0xb2bd0b28L, 0x2bb45a92L, 0x5cb36a04L, 0xc2d7ffa7L, 0xb5d0cf31L,
  1760.   0x2cd99e8bL, 0x5bdeae1dL, 0x9b64c2b0L, 0xec63f226L, 0x756aa39cL,
  1761.   0x026d930aL, 0x9c0906a9L, 0xeb0e363fL, 0x72076785L, 0x05005713L,
  1762.   0x95bf4a82L, 0xe2b87a14L, 0x7bb12baeL, 0x0cb61b38L, 0x92d28e9bL,
  1763.   0xe5d5be0dL, 0x7cdcefb7L, 0x0bdbdf21L, 0x86d3d2d4L, 0xf1d4e242L,
  1764.   0x68ddb3f8L, 0x1fda836eL, 0x81be16cdL, 0xf6b9265bL, 0x6fb077e1L,
  1765.   0x18b74777L, 0x88085ae6L, 0xff0f6a70L, 0x66063bcaL, 0x11010b5cL,
  1766.   0x8f659effL, 0xf862ae69L, 0x616bffd3L, 0x166ccf45L, 0xa00ae278L,
  1767.   0xd70dd2eeL, 0x4e048354L, 0x3903b3c2L, 0xa7672661L, 0xd06016f7L,
  1768.   0x4969474dL, 0x3e6e77dbL, 0xaed16a4aL, 0xd9d65adcL, 0x40df0b66L,
  1769.   0x37d83bf0L, 0xa9bcae53L, 0xdebb9ec5L, 0x47b2cf7fL, 0x30b5ffe9L,
  1770.   0xbdbdf21cL, 0xcabac28aL, 0x53b39330L, 0x24b4a3a6L, 0xbad03605L,
  1771.   0xcdd70693L, 0x54de5729L, 0x23d967bfL, 0xb3667a2eL, 0xc4614ab8L,
  1772.   0x5d681b02L, 0x2a6f2b94L, 0xb40bbe37L, 0xc30c8ea1L, 0x5a05df1bL,
  1773.   0x2d02ef8dL
  1774. };
  1775.  
  1776. unsigned long CRC32checksum(const char *seq, const long seqlen, unsigned long *checktotal)
  1777. /*CRC32checksum: modified from CRC-32 algorithm found in ZIP compression source */
  1778. {
  1779.   register unsigned long c = 0xffffffffL;
  1780.   register long n = seqlen;
  1781.     register char *sp = (char*) seq;
  1782.  
  1783.     while (n--) {
  1784.         c = crctab[((int)c ^ (to_upper(*sp))) & 0xff] ^ (c >> 8);
  1785.         sp++;
  1786.         }
  1787.   c= c ^ 0xffffffffL;
  1788.   *checktotal += c;
  1789.   return c;
  1790. }
  1791.  
  1792.  
  1793.  
  1794.  
  1795. short getseqtype( const char *seq, const long seqlen)
  1796. { /* return sequence kind: kDNA, kRNA, kProtein, kOtherSeq, ??? */
  1797.   char  c;
  1798.     register char *sp = (char*) seq;
  1799.   short i, maxtest;
  1800.   short na = 0, aa = 0, po = 0, nt = 0, nu = 0, ns = 0, no = 0;
  1801.  
  1802.   maxtest = (short) min(300, seqlen);
  1803.   for (i = 0; i < maxtest; i++) {
  1804. #if 1
  1805.         c =  *sp;  sp++;
  1806.     if (isprotonly(c)) po++;
  1807.     else if (isprimenuc(c)) {
  1808.       na++;
  1809.       if (isdnanuc(c)) nt++;
  1810.       else if (isrnanuc(c)) nu++;
  1811.       }
  1812.     else if (isamino(c)) aa++;
  1813.     else if (isseqsym(c)) ns++;
  1814.     else if (isalpha(c)) no++;
  1815. #else
  1816.     if (strchr(protonly, c)) po++;
  1817.     else if (strchr(primenuc,c)) {
  1818.       na++;
  1819.       if (c == 'T') nt++;
  1820.       else if (c == 'U') nu++;
  1821.       }
  1822.     else if (strchr(aminos,c)) aa++;
  1823.     else if (strchr(seqsymbols,c)) ns++;
  1824.     else if (isalpha(c)) no++;
  1825. #endif
  1826.     }
  1827.  
  1828.   if ((no > 0) || (po+aa+na == 0)) return kOtherSeq;
  1829.   /* ?? test for probability of kOtherSeq ?, e.g.,
  1830.   else if (po+aa+na / maxtest < 0.70) return kOtherSeq;
  1831.   */
  1832.   else if (po > 0) return kAmino;
  1833.   else if (aa == 0) {
  1834.     if (nu > nt) return kRNA;
  1835.     else return kDNA;
  1836.     }
  1837.   else if (na > aa) return kNucleic;
  1838.   else return kAmino;
  1839. } /* getseqtype */
  1840.  
  1841.  
  1842. void GetSeqStats( const char* seq, const long seqlen, 
  1843.                             short* seqtype, unsigned long* checksum, long* basecount)
  1844. {
  1845.     register unsigned long chk = 0xffffffffL;
  1846.   register long n = seqlen;
  1847.     register char c, *seqp;
  1848.     short         na = 0, aa= 0, no= 0, ns= 0, nt= 0, po= 0, nu= 0;
  1849.     short         seqt = kOtherSeq;
  1850.     unsigned long checks= 0;
  1851.     long     basec= 0;
  1852.     
  1853.     if (seq) {
  1854.         na= aa= no= ns= nt= po= nu= 0;
  1855.         seqp= (char*) seq;
  1856.         while (n--) {
  1857.             c= to_upper(*seqp);
  1858.             seqp++;
  1859.             if (c > ' '  && (c < '0' || c > '9')) { 
  1860.                 chk = crctab[((int)chk ^ c) & 0xff] ^ (chk >> 8);
  1861.             
  1862.                 basec++;
  1863. #if 1
  1864.             if (isprotonly(c)) po++;
  1865.             else if (isprimenuc(c)) {
  1866.               na++;
  1867.               if (isdnanuc(c)) nt++;
  1868.               else if (isrnanuc(c)) nu++;
  1869.               }
  1870.             else if (isamino(c)) aa++;
  1871.             else if (isseqsym(c)) ns++;
  1872.             else if (isalpha(c)) no++;
  1873. #else                
  1874.                 if (strchr( protonly, c)) po++;
  1875.                 else if (strchr(primenuc, c)) {
  1876.                     na++;
  1877.                     if (c == 'T') nt++;
  1878.                     else if (c == 'U') nu++;
  1879.                     }
  1880.                 else if (strchr(aminos, c))  aa++;  
  1881.                 else if (strchr(seqsymbols, c)) ns++;
  1882.                 else if (isalpha(c)) no++;
  1883. #endif
  1884.                 }
  1885.             }
  1886.       /* checks= chk ^ 0xffffffffL; */
  1887.         if (no > 0 || po+aa+na == 0) seqt = kOtherSeq;
  1888.         else if (po > 0) seqt = kAmino;
  1889.         else if (aa == 0) {
  1890.             if (nu > nt) seqt = kRNA;
  1891.             else seqt = kDNA;
  1892.             }
  1893.         else if (na > aa) seqt = kNucleic;
  1894.         else seqt = kAmino;
  1895.         }
  1896.     if (seqtype) *seqtype= seqt;
  1897.     if (checksum) *checksum= chk ^ 0xffffffffL;
  1898.     if (basecount) *basecount= basec;
  1899. }
  1900.  
  1901.  
  1902. char* compressSeq( const char gapc, const char *seq, const long seqlen, long *newlen)
  1903. {
  1904.   register char *a, *b;
  1905.   register long i;
  1906.   char  *newseq;
  1907.  
  1908.   *newlen= 0;
  1909.   if (!seq) return NULL;
  1910.   newseq = (char*) MALLOC(seqlen+1);
  1911.   if (!newseq) return NULL;
  1912.   for (a= (char*)seq, b=newseq, i=0; *a!=0; a++)
  1913.     if (*a != gapc) {
  1914.       *b++= *a;
  1915.       i++;
  1916.       }
  1917.   *b= '\0';
  1918.     REALLOC(newseq, i+1);
  1919.   *newlen= i;
  1920.   return newseq;
  1921. }
  1922.  
  1923.  
  1924.  
  1925. /***
  1926. char *rtfhead = "{\\rtf1\\defformat\\mac\\deff2 \
  1927. {\\fonttbl\
  1928.   {\\f1\\fmodern Courier;}{\\f2\\fmodern Monaco;}\
  1929.   {\\f3\\fswiss Helvetica;}{\\f4\\fswiss Geneva;}\
  1930.   {\\f5\\froman Times;}{\\f6\\froman Palatino;}\
  1931.   {\\f7\\froman New Century Schlbk;}{\\f8\\ftech Symbol;}}\
  1932. {\\stylesheet\
  1933.   {\\s1 \\f5\\fs20 \\sbasedon0\\snext1 name;}\
  1934.   {\\s2 \\f3\\fs20 \\sbasedon0\\snext2 num;}\
  1935.   {\\s3 \\f1\\f21 \\sbasedon0\\snext3 seq;}}";
  1936.  
  1937. char *rtftail = "}";
  1938. ****/
  1939.  
  1940.  
  1941.  
  1942.  
  1943.  
  1944.  
  1945. short writeSeq(FILE *outf, const char *seq, const long seqlen,
  1946.                 const short outform, const char *seqid)
  1947. {
  1948.     gFile = outf; gLinestart = 0;
  1949.     return writeSeqCall( readWriteFromFile, seq, seqlen, outform, seqid);
  1950. }
  1951.  
  1952.  
  1953. short writeSeqCall(
  1954.                                 ReadWriteProc writer, 
  1955.                                 const char *seq, const long seqlen,
  1956.                 const short outform, const char *seqid)
  1957. /* dump sequence to standard output */
  1958. {
  1959.   const short kSpaceAll = -9;
  1960. #define kMaxseqwidth  250
  1961.  
  1962.   boolean baseonlynum= false; /* nocountsymbols -- only count true bases, not "-" */
  1963.   short  numline = 0; /* only true if we are writing seqq number line (for interleave) */
  1964.   boolean numright = false, numleft = false;
  1965.   boolean nameright = false, nameleft = false, dumb = true;
  1966.   short   namewidth = 8, numwidth = 8;
  1967.   short   spacer = 0, width  = 50, tab = 0;
  1968.   /* new parameters: width, spacer, those above... */
  1969.     boolean dofreeseq= false;
  1970.     
  1971.   short linesout = 0, seqtype = kNucleic;
  1972.   long  i, j, l, l1, ibase;
  1973.   char  idword[31], endstr[30];
  1974.   char  seqnamestore[128], *seqname = seqnamestore;
  1975.   char  s[kMaxseqwidth], *cp;
  1976.   char  nameform[30], numform[30], nocountsymbols[30];
  1977.   unsigned long checksum = 0, checktotal = 0;
  1978.     char    outbuf[512];
  1979.  
  1980. #define FPUTC(c)    {outbuf[0]=c; outbuf[1]=0; (void)(*writer)(outbuf, 255, kRSFile_Write);}
  1981. #define FPUTS(s)    (void)(*writer)(outbuf, 255L, kRSFile_Write)
  1982.  
  1983.     /*fprintf(stderr,"writeSeqCall &writer= %ld  &gFile= %ld\n", writer, gFile);*/
  1984.  
  1985.   gPretty.atseq++;
  1986.   /* skipwhitespace(seqid); //???? */
  1987.   strncpy( seqnamestore, seqid, sizeof(seqnamestore));
  1988.   seqname[sizeof(seqnamestore)-1] = 0;
  1989.  
  1990.   sscanf( seqname, "%30s", idword);
  1991.   sprintf(numform, "%d", seqlen);
  1992.   numwidth= strlen(numform)+1;
  1993.   nameform[0]= '\0';
  1994.  
  1995.   if (strstr(seqname,"checksum") != NULL) {
  1996.     cp = strstr(seqname,"bases");
  1997.     if (cp!=NULL) {
  1998.       for ( ; (cp!=seqname) && (*cp!=','); cp--) ;
  1999.       if (cp!=seqname) *cp=0;
  2000.       }
  2001.     }
  2002.  
  2003.   strcpy( endstr,"");
  2004.   l1 = 0;
  2005.  
  2006.   if (outform == kGCG || outform == kMSF) {
  2007.       /* ! translate quasi-standard '-' indel to '.' that gcg uses */
  2008.       /* must do before checksum... */
  2009.       if (memchr(seq,'-',seqlen)) {
  2010.           char *newseq, *nc, *sc;
  2011.           newseq = (char*) MALLOC(seqlen+1);
  2012.             for (i=0, sc= (char*)seq, nc= newseq; i<=seqlen; i++, sc++, nc++) {
  2013.                 if (*sc == '-') *nc= '.'; 
  2014.                 else *nc= *sc;
  2015.                 }
  2016.             seq= (const char*) newseq;
  2017.             dofreeseq= true;
  2018.           }
  2019.     checksum = GCGchecksum(seq, seqlen, &checktotal);
  2020.     }
  2021.   else
  2022.     checksum = seqchecksum(seq, seqlen, &checktotal);
  2023.  
  2024.   switch (outform) {
  2025.  
  2026.     case kPlain:
  2027.     case kUnknown:    /* no header, just sequence */
  2028.       strcpy(endstr,"\n"); /* end w/ extra blank line */
  2029.       break;
  2030.  
  2031.     case kOlsen:  /* Olsen seqq. editor takes plain nucs OR Genbank  */
  2032.     case kGenBank:
  2033.       sprintf(outbuf,"LOCUS       %s       %d bp\n", idword, seqlen);  FPUTS(outbuf);
  2034.       sprintf(outbuf,"DEFINITION  %s  %d bases  %X checksum\n",  
  2035.                         seqname, seqlen, checksum); FPUTS(outbuf);
  2036.    /* sprintf(outbuf,"ACCESSION   %s\n", accnum); FPUTS(outbuf); */
  2037.       sprintf(outbuf,"ORIGIN      \n" ); FPUTS(outbuf);
  2038.       spacer = 11;
  2039.       numleft = true;
  2040.       numwidth = 8;  /* dgg. 1Feb93, patch for GDE fail to read short numwidth */
  2041.       strcpy(endstr, "\n//");
  2042.       linesout += 4;
  2043.       break;
  2044.  
  2045.     case kPIR:
  2046.       /* somewhat like genbank... \\\*/
  2047.       /* sprintf(outbuf,"\\\\\\\n"); FPUTS(outbuf); << only at top of file, not each entry... */
  2048.       sprintf(outbuf,"ENTRY           %s \n", idword); FPUTS(outbuf);
  2049.       sprintf(outbuf,"TITLE           %s  %d bases  %X checksum\n", 
  2050.                     seqname, seqlen, checksum); FPUTS(outbuf);
  2051.    /* sprintf(outbuf,"ACCESSION       %s\n", accnum); FPUTS(outbuf); */
  2052.       sprintf(outbuf,"SEQUENCE        \n" ); FPUTS(outbuf);
  2053.       numwidth = 7;
  2054.       width= 30;
  2055.       spacer = kSpaceAll;
  2056.       numleft = true;
  2057.       strcpy(endstr, "\n///");
  2058.       /* run a top number line for PIR */
  2059.       for (j=0; j<numwidth; j++) FPUTC(' ');
  2060.       for (j= 5; j<=width; j += 5) { sprintf(outbuf,"%10d", j ); FPUTS(outbuf); }
  2061.       FPUTC('\n');
  2062.       linesout += 5;
  2063.       break;
  2064.  
  2065.     case kNBRF:
  2066.       if (getseqtype(seq, seqlen) == kAmino)
  2067.         { sprintf(outbuf,">P1;%s\n", idword );  FPUTS(outbuf);}
  2068.       else
  2069.         { sprintf(outbuf,">DL;%s\n", idword );  FPUTS(outbuf);}
  2070.       sprintf(outbuf,"%s  %d bases  %X checksum\n", 
  2071.                 seqname, seqlen, checksum); FPUTS(outbuf);
  2072.       spacer = 11;
  2073.       strcpy(endstr,"*\n");
  2074.       linesout += 3;
  2075.       break;
  2076.  
  2077.     case kEMBL:
  2078.       sprintf(outbuf,"ID   %s\n", idword ); FPUTS(outbuf);
  2079.   /*  sprintf(outbuf,"AC   %s\n", accnum ); FPUTS(outbuf); */
  2080.       sprintf(outbuf,"DE   %s  %d bases %X checksum\n", 
  2081.                         seqname, seqlen, checksum); FPUTS(outbuf);
  2082.       sprintf(outbuf,"SQ             %d BP\n", seqlen ); FPUTS(outbuf);
  2083.       strcpy(endstr, "\n//"); /* 11Oct90: bug fix*/
  2084.       tab = 4;     /** added 31jan91 */
  2085.       spacer = 11; /** added 31jan91 */
  2086.       width = 60;
  2087.       linesout += 4;
  2088.       break;
  2089.  
  2090.     case kGCG:
  2091.       sprintf(outbuf,"%s\n", seqname ); FPUTS(outbuf);
  2092.    /* sprintf(outbuf,"ACCESSION   %s\n", accnum );  FPUTS(outbuf);*/
  2093.       sprintf(outbuf,"    %s  Length: %d  (today)  Check: %d  ..\n", 
  2094.                 idword, seqlen, checksum ); FPUTS(outbuf);
  2095.       spacer = 11;
  2096.       numleft = true;
  2097.       strcpy(endstr, "\n");  /* this is insurance to help prevent misreads at eof */
  2098.       linesout += 3;
  2099.       break;
  2100.  
  2101.     case kStrider: /* ?? map ?*/
  2102.       sprintf(outbuf,"; ### from DNA Strider ;-)\n" ); FPUTS(outbuf);
  2103.       sprintf(outbuf,"; DNA sequence  %s  %d bases  %X checksum\n;\n", 
  2104.                 seqname, seqlen, checksum); FPUTS(outbuf);
  2105.       strcpy(endstr, "\n//");
  2106.       linesout += 3;
  2107.       break;
  2108.  
  2109.     case kFitch:
  2110.       sprintf(outbuf,"%s  %d bases  %X checksum\n", 
  2111.                 seqname, seqlen, checksum); FPUTS(outbuf);
  2112.       spacer = 4;
  2113.       width = 60;
  2114.       linesout += 1;
  2115.       break;
  2116.  
  2117.     case kPhylip2:
  2118.     case kPhylip4:
  2119.       /* this is version 3.2/3.4 -- simplest way to write
  2120.         version 3.3 is to write as version 3.2, then
  2121.         re-read file and interleave the species lines */
  2122.       if (strlen(idword)>10) idword[10] = 0;
  2123.       sprintf(outbuf,"%-10s  ", idword ); FPUTS(outbuf);
  2124.       l1  = -1;
  2125.       tab = 12;
  2126.       spacer = 11;
  2127.       break;
  2128.  
  2129.     case kASN1:
  2130.       seqtype= getseqtype(seq, seqlen);
  2131.       switch (seqtype) {
  2132.         case kDNA     : cp= "dna"; break;
  2133.         case kRNA     : cp= "rna"; break;
  2134.         case kNucleic : cp= "na"; break;
  2135.         case kAmino   : cp= "aa"; break;
  2136.         case kOtherSeq: cp= "not-set"; break;
  2137.         }
  2138.       sprintf(outbuf,"  seq {\n" ); FPUTS(outbuf);
  2139.       sprintf(outbuf,"    id { local id %d },\n", gPretty.atseq ); FPUTS(outbuf);
  2140.       sprintf(outbuf,"    descr { title \"%s\" },\n",  seqid ); FPUTS(outbuf);
  2141.       sprintf(outbuf,"    inst {\n", dumb ); FPUTS(outbuf);
  2142.       sprintf(outbuf,"      repr raw, mol %s, length %d, topology linear,\n", 
  2143.                 cp, seqlen ); FPUTS(outbuf);
  2144.       sprintf(outbuf,"      seq-data\n" ); FPUTS(outbuf);
  2145.       if (seqtype == kAmino)
  2146.         { sprintf(outbuf,"        iupacaa \""); FPUTS(outbuf);}
  2147.       else
  2148.         { sprintf(outbuf,"        iupacna \"" ); FPUTS(outbuf);}
  2149.       l1  = 17;
  2150.       spacer = 0;
  2151.       width  = 78;
  2152.       tab  = 0;
  2153.       strcpy(endstr,"\"\n      } } ,");
  2154.       linesout += 7;
  2155.       break;
  2156.  
  2157.     case kPAUP:
  2158.       nameleft= true;
  2159.       namewidth = 9;
  2160.       spacer = 21;
  2161.       width  = 100;
  2162.       tab  = 0; /* 1; */
  2163.       /* strcpy(endstr,";\nend;"); << this is end of all seqs.. */
  2164.       /* do a header comment line for paup */
  2165.       sprintf(outbuf,"[Name: %-16s  Len:%6d  Check: %8X]\n", 
  2166.                 idword, seqlen, checksum ); FPUTS(outbuf);
  2167.       linesout += 1;
  2168.       break;
  2169.  
  2170.     case kPretty:
  2171.       numline= gPretty.numline;
  2172.       baseonlynum= gPretty.baseonlynum;
  2173.       namewidth = gPretty.namewidth;
  2174.       numright = gPretty.numright;
  2175.       numleft = gPretty.numleft;
  2176.       nameright = gPretty.nameright;
  2177.       nameleft = gPretty.nameleft;
  2178.       spacer = gPretty.spacer + 1;
  2179.       width  = gPretty.seqwidth;
  2180.       tab  = gPretty.tab;
  2181.       /* also add rtf formatting w/ font, size, style */
  2182.       if (gPretty.nametop) {
  2183.         sprintf(outbuf,"Name: %-16s  Len:%6d  Check: %8X\n", 
  2184.                     idword, seqlen, checksum ); FPUTS(outbuf);
  2185.         linesout++;
  2186.         }
  2187.       break;
  2188.  
  2189.     case kMSF:
  2190.       sprintf(outbuf," Name: %-16s Len:%6d  Check:%5d  Weight:  1.00\n",
  2191.                      idword, seqlen, checksum ); FPUTS(outbuf);
  2192.       linesout++;
  2193.       nameleft= true;
  2194.       namewidth= 15; /* need MAX namewidth here... */
  2195.       sprintf(nameform, "%%+%ds ",namewidth);
  2196.       spacer = 11;
  2197.       width  = 50;
  2198.       tab = 0; /* 1; */
  2199.       break;
  2200.  
  2201.     case kIG:
  2202.       sprintf(outbuf,";%s  %d bases  %X checksum\n",  
  2203.                 seqname, seqlen, checksum ); FPUTS(outbuf);
  2204.       sprintf(outbuf,"%s\n",  idword ); FPUTS(outbuf);
  2205.       strcpy(endstr,"1"); /* == linear dna */
  2206.       linesout += 2;
  2207.       break;
  2208.  
  2209.     default :
  2210.     case kZuker: /* don't attempt Zuker's ftn format */
  2211.     case kPearson:
  2212.       sprintf(outbuf,">%s  %d bases  %X checksum\n",  
  2213.                 seqname, seqlen, checksum ); FPUTS(outbuf);
  2214.       linesout += 1;
  2215.       break;
  2216.     }
  2217.  
  2218.   if (*nameform==0) sprintf(nameform, "%%%d.%ds ",namewidth,namewidth);
  2219.   if (numline) sprintf(numform, "%%%ds ",numwidth);
  2220.   else sprintf(numform, "%%%dd ",numwidth);
  2221.   strcpy( nocountsymbols, kNocountsymbols);
  2222.   if (baseonlynum) {
  2223.     if (strchr(nocountsymbols,gPretty.gapchar)==NULL) {
  2224.       strcat(nocountsymbols," ");
  2225.       nocountsymbols[strlen(nocountsymbols)-1]= gPretty.gapchar;
  2226.       }
  2227.     if (gPretty.domatch && (cp=strchr(nocountsymbols,gPretty.matchchar))!=NULL) {
  2228.       *cp= ' ';
  2229.       }
  2230.     }
  2231.  
  2232.   if (numline) {
  2233.    *idword= 0;
  2234.    }
  2235.  
  2236.   width = min(width,kMaxseqwidth);
  2237.   for (i=0, l=0, ibase = 1; i < seqlen; ) {
  2238.  
  2239.     if (l1 < 0) l1 = 0;
  2240.     else if (l1 == 0) {
  2241.       if (nameleft) { sprintf(outbuf, nameform,  idword ); FPUTS(outbuf);}
  2242.       if (numleft) { if (numline) { sprintf(outbuf, numform,  "" ); FPUTS(outbuf);}
  2243.                     else { sprintf(outbuf, numform, ibase );  FPUTS(outbuf);}
  2244.                                         }
  2245.       for (j=0; j<tab; j++) FPUTC(' ');
  2246.       }
  2247.  
  2248.     l1++;                 /* don't count spaces for width*/
  2249.     if (numline) {
  2250.       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1)) {
  2251.         if (numline==1) FPUTC(' ');
  2252.         s[l++] = ' ';
  2253.         }
  2254.       if (l1 % 10 == 1 || l1 == width) {
  2255.         if (numline==1) { sprintf(outbuf,"%-9d ", i+1 ); FPUTS(outbuf);}
  2256.         s[l++]= '|'; /* == put a number here */
  2257.         }
  2258.       else s[l++]= ' ';
  2259.       i++;
  2260.       }
  2261.  
  2262.     else {
  2263.       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
  2264.         s[l++] = ' ';
  2265.       if (!baseonlynum) ibase++;
  2266.       else if (0==strchr(nocountsymbols,seq[i])) ibase++;
  2267.       s[l++] = seq[i++];
  2268.       }
  2269.  
  2270.     if (l1 == width || i == seqlen) {
  2271.       if (outform==kPretty) for ( ; l1<width; l1++) {
  2272.         if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
  2273.           s[l++] = ' ';
  2274.         s[l++]=' '; /* pad w/ blanks */
  2275.         }
  2276.       s[l] = '\0';
  2277.       l = 0; l1 = 0;
  2278.  
  2279.       if (numline) {
  2280.         if (numline==2) { sprintf(outbuf,"%s", s );  FPUTS(outbuf);} /* finish numberline ! and | */
  2281.         }
  2282.       else {
  2283.         if (i == seqlen) { sprintf(outbuf,"%s%s", s,endstr ); FPUTS(outbuf); }
  2284.         else { sprintf(outbuf,"%s", s ); FPUTS(outbuf);}
  2285.         if (numright || nameright) FPUTC(' ');
  2286.         if (numright)  { sprintf(outbuf,numform,  ibase-1 ); FPUTS(outbuf);}
  2287.         if (nameright) { sprintf(outbuf, nameform, idword ); FPUTS(outbuf);}
  2288.         }
  2289.       FPUTC('\n');
  2290.       linesout++;
  2291.       }
  2292.     }
  2293.   if (dofreeseq) FREE((char*)seq);
  2294.   return linesout;
  2295.     
  2296. #undef FPRINTF 
  2297. #undef FPUTC 
  2298. }  /*writeSeq*/
  2299.  
  2300.  
  2301.  
  2302. /* End file: ureadseq.c */
  2303.